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Abstract 


In a direct interplanetary transfer, the spacecraft moves from a parking orbit of the departure planet to a parking orbit of the arrival 
planet. The transfer trajectory must be designed such that the specified arrival parking orbit conditions are achieved. For a fixed depar- 
ture epoch and flight duration, there are four distinct transfer trajectory design options in a direct transfer. The conventional patched 
conic method, the most widely used analytical trajectory design method, does not identify these design options. An iterative patched 
conic method that identifies these distinct design options is developed and presented in this paper. This method involves two iterative 
processes: (i) iteration on the hyperbolic orbit characteristics using an analytical tuning strategy to achieve the hyperbolic excess velocity 
vector at the patch point, (ii) iteration on the patch points at the sphere of influence. The performance of the proposed method is com- 
pared with the conventional and V-infinity tuned patched conic methods. A design analysis tool, based on the proposed method, is devel- 


oped and tested in various orbiter mission scenarios. 
© 2017 COSPAR. Published by Elsevier Ltd. All rights reserved. 


Keywords: Orbiter mission; Patched conic; Sphere of influence; V-infinity 


1. Introduction 


In a direct interplanetary transfer, the spacecraft experi- 
ences the gravitational pull of multiple celestial bodies. The 
unknown transfer trajectory characteristics are obtained by 
numerically solving the n-body problem using Cowell’s 
method. However, the numerical method requires an initial 
guess. The analytical methods generate designs quickly 
which can be used as initial guess for numerical refinement. 
Walter Hohmann (1960) introduced an analytical method 
which is essentially a heliocentric transfer between the orbits 
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of the target planets assuming that the planetary orbits are 
circular and coplanar. The interplanetary trajectory design 
using Hohmann method generates only the magnitude of 
the velocity impulses and does not generate the direction. 
Also, the actual positions and gravity fields of the planets 
are not considered. The conventional patched conic method 
overcomes these difficulties and provides improved design. 
The preliminary mission design and analysis, in general, is 
carried out using the conventional patched conic method 
(Heaton et al., 2002; Englander et al., 2012; Spreen et al., 
2011; Campagnola et al., 2014). The linked conic method 
further improves the design (Cornelisse, 1978). These meth- 
ods are briefly discussed below. 


1.1. Conventional patched conic method 


The conventional patched conic method considers the 
target planets as points in space and connects them through 
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Nomenclature 

a semi-major axis, km 

e eccentricity 

i inclination, deg 

R heliocentric position vector, km 

T planetocentric position vector, km 

ta SOI duration for the arrival planet, days 

tp SOI duration for the departure planet, days 
trp flight duration, days 

Vv heliocentric velocity vector, km/s 

v planetocentric velocity vector, km/s 

Y angle between the asymptotic velocity vector 


and the propagated velocity vector, deg 


Q right ascension of ascending node, deg 

Olas right ascension of the asymptotic velocity vec- 
tor, deg 

O00 declination of the asymptotic velocity vector, 
deg 

€ threshold value 

Dx asymptotic true anomaly, deg 

m gravitational constant, km?/s* 


v true anomaly, deg 
o argument of periapsis, deg 
Subscripts 


A arrival planet 

D departure planet 
P periapsis 

T transfer trajectory 
oe) hyperbolic orbit 


Abbreviations 
argument of periapsis 


APO arrival parking orbit 
CAA closest approach altitude 
DPO departure parking orbit 


PCFM patched conic force model 

POI planetary orbit insertion 

RAAN right ascension of ascending node 
SOI sphere of influence 

trajectory correction maneuver 
TPI trans-planetary injection 


a heliocentric Lambert conic (Battin, 1987). The gravita- 
tional forces of the target planets are considered and super- 
imposed within the respective sphere of influences (SOI). 
Accordingly, the transfer trajectory is subdivided into three 
two-body trajectories: (i) a planetocentric departure hyper- 
bola till the SOI of the departure planet, (11) a heliocentric 
ellipse from the SOI of the departure planet till the SOI of 
the arrival planet, and (iii) a planetocentric arrival hyper- 
bola within the SOI of the arrival planet. Because the helio- 
centric velocities in the transfer trajectory at the departure 
and SOI differ by less than 1 m/s, the excess velocity vector 
is computed as the difference between the Lambert velocity 
vector and the planet velocity vector at the departure 
epoch. The planetocentric hyperbolic orbit is obtained 
from the excess velocity vector. In this technique, there is 
discontinuity in the velocity vectors at the SOI because of 
the superimposition of gravity of the planet within the 
SOI. A hybrid patched conic method is discussed in the 
context of lunar transfer (Escobal, 1968). 


1.2. Linked conic method 


In this method, the planetocentric and heliocentric tra- 
jectories are patched at the SOI using different strategies 
to ensure continuity. Clarke et al. (1966) used the differen- 
tial correction method to modify the departure excess 
velocity vectors such that the desired arrival B-plane coor- 
dinates are achieved. Cornelisse (1978) synchronized the 
velocity at the end of the geocentric phase with the velocity 
at the beginning of the heliocentric phase using Newton- 


Raphson method. Ramanan and Adimurthy (2005) used 
an analytical tuning strategy in the context of lunar trans- 
fers to patch up the velocity discontinuities at the SOI 
boundary. The departure hyperbolic orbit characteristics 
are modified such that the asymptotic excess velocity vector 
is achieved at the SOI. Parvathi and Ramanan (2016) used 
this analytical tuning strategy to tune the hyperbolic orbit 
characteristics in the interplanetary context and termed the 
method as V-infinity tuned patched conic method. 


1.3. Motivation for the current work 


Both the conventional and V-infinity tuned patched 
conic methods cannot generate the four distinct transfer 
trajectory design options that exist for a direct orbiter mis- 
sion opportunity. Sergeyevsky et al. (1983) applied the 
pseudostate technique that was introduced by Wilson 
(1970), for direct interplanetary transfers. But, rectilinear 
hyperbolas were used for the planetocentric phases. 
Parvathi and Ramanan (2016) presented an iterative pseu- 
dostate method by assuming non-rectilinear hyperbolas 
that can identify the distinct design options. In this paper, 
an iterative analytical method based on patched conic tech- 
nique is introduced considering the complexity involved in 
the implementation of the iterative pseudostate method. 
The proposed method generates the distinct design options 
while retaining the simplicity of the patched conic approx- 
imations. Further, in the absence of additional information 
on the arrival geometry, the numerical refinement of the 
design obtained using the conventional and V-infinity 
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tuned patched conic methods can converge to any one of 
the four distinct design options. The additional informa- 
tion on arrival angles, Right Ascension of Ascending Node 
(RAAN) and Argument of Periapsis (AoP), can be 
obtained from the iterative patched conic method. 

The paper is organized as follows. The method to gener- 
ate the orbit characteristics from the excess velocity vector 
is described in Section 2. The analytical tuning strategy, 
used to modify the hyperbolic orbit characteristics that 
achieve the excess velocity vector at the patch point on 
SOI, is described in Section 3. The iterative patched conic 
method that iterates on the patch point is discussed in 
detail in Section 4. The numerical results for various orbi- 
ter mission scenarios are analyzed in Section 5. 


2. Determination of hyperbolic/parking orbit characteristics 


The shape of the planets are assumed to be spherical. 
The hyperbolic orbit characteristics are obtained from the 
excess velocity vector by assuming a periapsis altitude 
and inclination. To minimize the trans-planetary injection 
(TPI) and planetary orbit insertion (POI) velocity impulses, 
it is essential that (i) the parking orbit plane is coplanar 
with the hyperbolic orbit plane, (ii) the location of TPI/ 
POI is at the periapsis of the parking orbit and (iii) the 
argument of periapsis (AoP) of the hyperbolic trajectory 
is the same as the argument of periapsis of the elliptical 
parking orbit or the argument of latitude of the circular 
parking orbit. Hence, the parking orbit characteristics 
except the apoapsis altitude are same as the hyperbolic 
orbit characteristics. The aforementioned conditions lead 
to tangential and horizontal injection of the velocity 
impulse. 

The geometry of the departure hyperbolic orbit/depar- 
ture parking orbit (DPO) is given in Fig. |. The departure 
asymptotic excess velocity vector (Vp) is represented by 
spherical coordinates (Vp, %.0p,0x0p) Which represents 
the magnitude, right ascension and declination of Vp 
respectively. The departure hyperbolic orbit characteristics 


Equator 


Y (x axis) 
Plane of parking orbit 
and Transfer orbit 


Nodal line 


Fig. 1. Geometry of departure hyperbolic orbit/DPO (Parvathi and 
Ramanan, 2016). 


are obtained from the departure asymptotic excess velocity 
vector using the following relations. The semi-major axis of 
the departure hyperbolic orbit is given by: 


Axop = —Hp/Veop (1) 


where lp is the gravitational constant of the departure pla- 
net. The related eccentricity is given by: 


Coop = Lae (TP Vocp/Hp) (2) 


where rp_,, is the magnitude of position vector at the peri- 
apsis of the departure hyperbolic orbit. 

The orientation of the departure hyperbolic orbit/DPO 
plane with respect to the planetocentric equatorial plane 
is given by two angles, RAAN (Q,.p) and inclination 
(ixp). The locations of the periapsis and the asymptotic 
radius vector are given by wp and 0..p respectively. In a 
departure hyperbola, the asymptotic radius vector and 
the asymptotic excess velocity vector are in the same direc- 
tion. Therefore, the direction of asymptotic excess velocity 
vector (V.op) is represented by the direction of the asymp- 
totic radius vector (fp) in Fig. 1. The departure angles, 
RAAN and AoP of the departure hyperbolic orbit are 
derived using spherical trigonometric relations (Tolson, 
1963; Ramanan, 2002). 


0...) = tito, tay (3) 


SIN(Moop + Oop) = SiN b.0p/ SIN in, (4) 


SIN(A0op — 


where 0p is given by, 
Osan = 008 1(—1/e,). (5) 


The value of inclination (i,.p) must be greater than the 
declination of the departure excess velocity vector to ensure 
the coplanarity of the parking orbit and the hyperbolic tra- 
jectory planes. 


Psy OS (6) 


The Eqs. (3) and (4) have two solutions respectively, viz. 


Qeop|1 = Mop — Sin ' (tan Soop / taN top) (7) 
Qoop | = %oop — 180 + sin! (tan 5,0, / tan irop) (8) 
Wvon|, = Sin "(SiN dcon/ Si ing, ) — Poon (9) 
Ooop|> = 180 — sin! (sin b.c,/ SiN icp) — Doop (10) 


These solutions correspond to the two possible geome- 
tries of departure hyperbolic orbit planes containing the 
excess velocity vector (Vp) while (i) ascending (ii) descend- 
ing. The Eqs. (3) and (4) have mathematical singularity 
when i,,,, = 0. However, in such a case, the declination of 
the excess velocity vector (6,,,) must also be zero. That 
means, i, = d.p, The solution can be obtained using this 
relation in the above equations. Note that for this case, 
there is only one geometry for the hyperbolic orbit plane. 

The two possible arrival hyperbolic orbit planes are 
obtained from the above equations by replacing the sub- 
script ‘D’ with ‘A’. The arrival asymptotic excess velocity 
vector (V..4) is represented using the spherical coordinates 
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Table 1 

Transfer trajectory design options (Parvathi and Ramanan, 2016). 
Transfer trajectory design options Departure Arrival 
Option 1-1 (Qoop 1; Poop 11) (Qo04|1; Poo, |1) 
Option 1-2 (Qson|1; @con 1) (Qo0,|2, Moo, |2) 
Option 2-1 (Qxp|2; Oop!) (Qo04|15 Poo, |1) 
Option 2-2 (Qn lo, Boop |>) (Qoo4 lo, Vv, |>) 


(VocA; %ccd; 9004) Which represents the magnitude, right 
ascension and declination of V4. In an arrival hyperbola, 
the asymptotic radius vector and the asymptotic excess 
velocity vector are in opposite directions. Hence, for the 
computation of the arrival hyperbolic orbit characteristics, 
the parameters (V0.4, %oc4; Ooo) are set as, 


sod = T+ Oso4 (11) 
Oocd = —Oo0A (12) 


Each of the departure hyperbolic orbit can be mapped 
to each of the arrival hyperbolic orbit, thus resulting in 
four transfer trajectory design options, given in Table 1. 


3. Analytical tuning strategy 


The departure hyperbolic orbit characteristics must be 
modified such that the asymptotic excess velocity vector 
(named as desired excess velocity vector) is achieved at 
the SOI upon Keplerian propagation. In Keplerian propa- 
gation, a state is obtained in one step based on two-body 
spherical force model. In this paper, the SOI is defined 
by the number of days a spacecraft spends within it. The 
initial estimate of the departure hyperbolic orbit character- 
istics is obtained from the departure asymptotic excess 
velocity vector using the procedure described in Section 2. 
The Keplerian propagation of the initial departure hyper- 
bolic orbit characteristics from the DPO periapsis for 
SOI duration, results in the propagated excess velocity vec- 
tor (Vap). This will not be the same as the desired excess 
velocity vector (Vp). The angle between ¥,.p and Vpp is 
represented by '¥ (cf. Fig. 2). In order to achieve the desired 
excess velocity vector at the SOI, the departure hyperbolic 


LV qo 


orbit characteristics are tuned using an analytical tuning 
strategy adapted from Ramanan and Adimurthy (2005). 

The philosophy of analytical tuning strategy in the inter- 
planetary context is described in Parvathi and Ramanan 
(2016). The propagated excess velocity vector is rotated 
by the angle (’) about the angular momentum vector keep- 
ing the DPO periapsis altitude fixed. The rotation causes a 
shift in the location of DPO periapsis. The magnitude of 
the desired excess velocity vector is achieved by changing 
the semi-major axis of the departure hyperbolic orbit. 
Because the DPO periapsis altitude is kept as the same, 
the eccentricity of the departure hyperbola changes. This 
process is repeated till the propagated excess velocity vec- 
tor matches the desired excess velocity vector at SOI. Note 
that the plane of the departure hyperbola does not change 
during this process. For completion, the tuning strategy for 
the departure phase is outlined as follows. 


(1) Fix the periapsis (rp_,,) and inclination (ip) of the 
departure hyperbolic orbit. 

(2) Compute the semi-major axis (a) and the related 
eccentricity (€.») of the departure hyperbolic orbit 
using Eqs. (1) and (2). 

(3) Compute the two values of RAAN (Q,,p) and AoP 
(Mp) from Eqs. (7)-(10), and choose one of the 
two departure hyperbolic geometries. 

(4) Propagate forward the departure hyperbolic orbit 
characteristics (Assp, €cx0D; hood; Qed; OooDs VooD = VP.) = 0) 
for a prefixed SOI duration in one step under the two- 
body Keplerian force model. Compute the position vector 
(Tap) and the propagated excess velocity vector (Vpp). 

(5) Compute € = |Vxp — Vap|. If € is less than a prefixed 
threshold value, the departure hyperbolic orbit char- 
acteristics are achieved, otherwise continue with the 
following steps. 

(6) Find the angle (‘¥) between the desired and propa- 
gated excess velocity vector viz., (Von) and (Vap). 

(7) The shift in AoP location is computed as: 


_-*”” Desired departure 
yet excess velocity vector 


Tuned 
Perigee 


*RAAN and Inclination do not change 


Fig. 2. Analytical tuning strategy at the departure end (Parvathi and Ramanan, 2016). 


S.P. Parvathi, R.V. Ramanan! Advances in Space Research 59 (2017) 1763-1774 1767 


(8) Compute the new location of periapsis of the depar- 
ture hyperbolic trajectory (a,,,) as: 


Op = Oxop + Aw (14) 


(9) Compute the new semi-major axis of the rotated 
departure hyperbolic trajectory from the equation: 


2 
dxop = “tp / (v0 = Ho) (15) 
lap 


and the related eccentricity from Eq. (2). 
(10) Repeat the steps (4)-(10) till the threshold value is 
achieved. 


For the arrival phase, the Eqs. (13)-(15) are used for 
analytical tuning by replacing the subscript ‘D’ by ‘A’. 
The steps followed are also same except step 4. The arrival 
hyperbolic orbit characteristics are propagated backward 
in one step for a prefixed arrival SOI duration under the 
two-body force model. 


4. Iterative patched conic method 


The iterative patched conic method identifies the distinct 
transfer trajectory design options available for an inter- 
planetary orbiter mission while maintaining the simplicity 
of the conventional patched conic method. The proposed 
method is applicable for direct transfer to all planets. For 
demonstration purposes, the three segments of a typical 
Earth to Mars transfer arising based on patched conic 
assumptions are given in Fig. 3. 

The working principle of the iterative patched conic 
method is as follows. Initially, the position vectors of the 
target planets, Earth and Mars, at the departure and arri- 
val epochs respectively, are connected using the heliocen- 
tric Lambert conic. The heliocentric velocities of the 
spacecraft in the transfer conic at the target point and at 
the SOI are assumed to be same. Therefore, the initial 
asymptotic excess velocity vector is calculated as the differ- 
ence between the heliocentric velocity vector in the transfer 


% \ 
; ? 
. A \ 
Heliocentric ellipse = 7 - ‘ ‘ 
. ~~ transfer trajectory \ B 
' ' Vay ‘ ; 
: >: ; ! 
SUN ! ; 
v7 . 
. \ , i 
Marscentric hyperbola ; ra ’ 
. a ‘ 
see 2 Z 
ars SOI ri 


Si es one ae? 


Fig. 3. Patched conic approximation for Earth-Mars direct interplanetary 
transfer. 


conic and the heliocentric velocity vector of the planet at 
the departure epoch. The initial hyperbolic orbit character- 
istics are computed from the excess velocity vector and they 
are modified using the analytical tuning strategy (cf. Sec- 
tion 3) to achieve the asymptotic excess velocity vector at 
the SOI. The planetocentric position vectors at the SOI 
are obtained by Keplerian propagation of the modified 
hyperbolic orbit characteristics. These locations are known 
as the patch points. The planetocentric position vectors of 
the patch points are transformed to the heliocentric frame 
and the patch points are connected using a heliocentric 
Lambert conic. The flight duration for the Lambert conic 
excludes the SOI durations. The excess velocity vector is 
updated and the hyperbolic orbit characteristics are modi- 
fied to obtain the updated excess velocity vector at the SOI. 
There are two loops involved, (i) the inner loop modifies 
the hyperbolic orbit characteristics to achieve the excess 
velocity vector at the SOI, (ii) the outer loop iterates on 
the patch points. For better understanding, the iterative 
patched conic algorithm is explained as follows. 


(1) For a departure epoch and flight duration (tgp), fix 
the periapsis altitude and inclination of the departure 
and arrival hyperbolic orbits. 

(2) Fix the SOI durations for the departure and arrival 
planets as tp and ta respectively. 

(3) Obtain the heliocentric states of the departure planet 
on the departure epoch (Rp, Vp) and the arrival pla- 
net on the arrival epoch (Ra, Va) from JPL ephemeris 
DE405. These heliocentric states are in Earth Mean 
Equator and Equinox of J2000 frame. 

(4) Connect the position vectors (Rp and Ra) through 
the Lambert conic for the flight duration (trp). Let 
(Vr, and V;,) be the corresponding heliocentric 
velocity vectors in the transfer trajectory. 

(5) Compute the planetocentric excess velocity vectors of 
the spacecraft at the departure and arrival ends, Vx, 
and Vx, respectively: 


=Vr, - Vp (16) 
Vac, = Wr, — Va (17) 


These planetocentric excess velocity vectors are in 
Earth Mean Equator and Equinox of J2000 frame. 
The arrival excess velocity vector is transformed to 
Mars Mean Equator and Equinox of J2000 (v’..,) 
(Archinal et al., 2009). 

(6) Compute the hyperbolic orbit characteristics 
(x05 C05 foo, Qo Woo} Voo = Ve,, = 0) of the departure 
and arrival hyperbolic orbits from the excess velocity 
vectors using the procedure given in Section 2. 


The steps above constitute the conventional patched 
conic method and result in four design options which are 
not distinct. One of the four design options is chosen for 
further modification. 
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(7) Propagate forward the departure hyperbolic orbit 
characteristics for SOI duration (tp). Similarly, prop- 
agate backward the arrival hyperbolic orbit charac- 
teristics for SOI duration (ta) from APO periapsis. 
These Keplerian propagations are under the two- 
body force model Keplerian. Let the propagated 
planetocentric velocity vectors at the departure and 
arrival SOIs be, ¥p,, and v’/,, respectively. 

(8) Use the analytical tuning strategy, described in Sec- 
tion 3, to match the desired and the propagated 
excess velocity vectors at the SOI: V., and v}, at 
departure, and v’,,, and v/,, on arrival respectively. 

(9) Obtain the planetocentric position vectors at SOI, fr, 
at departure and r’,, on arrival, by propagating the 
respective tuned hyperbolic orbit characteristics 
under two-body Keplerian force model. These are 
the patch points. The planetocentric position vector 
of the arrival patch point (1’,,) is transformed to 
Earth equator and Equinox of J2000 frame (f,,) 
because the tuned arrival hyperbolic orbit character- 
istics are in Mars centered coordinate frame. 


Step 7 involves the numerical solution of Kepler equa- 
tion. For a threshold value of 1E—15, each solution run 
takes about 4—5 iterations. The tuning strategy is termed 
as ‘analytical’ because it does not involve numerical inte- 
gration of the equations of motion. The steps 1-9 consti- 
tute the V-infinity tuned patched conic method. 


(10) Transform the planetocentric position vectors of the 
patch points at the departure and arrival SOI into 
heliocentric position vectors (RS! and R$), 


RE — RE +3, (18) 
RQ = Ra +Tha (19) 


where Rj is the heliocentric position vector of the 
departure planet after tp days from the departure 
epoch and Ri is the heliocentric position vector of 
the arrival planet before ta days from the arrival 
epoch. These position vectors are with respect to 
Earth Equator and Equinox of J2000 frame. 

(11) Connect the heliocentric patch point position vectors 
(R$! and R§!) using the heliocentric Lambert conic 
for the flight duration (trp — tp — ta). Let the corre- 
sponding heliocentric velocity vectors at the patch 
points be V7°! and V70! respectively. 

(12) The planetocentric velocity vectors of the spacecraft 
at the patch points, Vp and V4 respectively, are 
now obtained as: 


Vay = VSO! — Wp (20) 
Vay = VIO VE (21) 


where Vj? is the heliocentric velocity vector of the 
departure planet after tp days from the departure 


epoch and Vi’ is the heliocentric velocity vector of 
the arrival planet before ta days from the arrival 
epoch. These velocity vectors are with respect to 
Earth Equator and Equinox of J2000 frame. The arri- 
val excess velocity vector is transformed to Mars 
Equator and Equinox of J2000 frame. 

(13) Compute the hyperbolic orbit characteristics 
(4x5 C504 bo01 250; Boos Yoo = Vp,, = 0) of departure and 
arrival hyperbolic orbits from the excess velocity vec- 
tors (cf. Section 2). 

(14) Propagate the hyperbolic orbit characteristics analyt- 
ically under two-body Keplerian force model for the 
SOI durations. The hyperbolic orbit characteristics 
are synchronized using the analytical tuning strategy 
(cf. Section 3) to achieve the desired excess velocity 
vector at the SOI. 

(15) Obtain the planetocentric position vector of the 
departure and arrival patch points. The planetocen- 
tric position vector of the arrival patch point is trans- 
formed to Earth Equator and Equinox of J2000 
frame. 

(16) Transform the planetocentric position vector of the 
departure and arrival patch points to heliocentric 
position vectors. The difference in magnitude between 
the heliocentric position vectors of two successive 
patch points are computed. 

(17) If the differences are less than the respective threshold 
values, then the transfer trajectory design is obtained. 
Otherwise the heliocentric position vectors of the 
patch points are updated with the new values and 
the steps 11-17 are repeated. 


At the end of these steps, the distinct transfer trajectory 
design for one of the options is obtained. We can apply 
these steps to obtain other three distinct designs as well. 


5. Results 
5.1. Performance comparison 


A FORTRAN code has been developed based on the 
proposed iterative patched conic method and used to gen- 
erate the distinct design options for a typical direct Earth 
to Mars orbiter mission. The earth parking orbit is 
assumed as 300 x 25,000 km and inclined at 75 degree with 
respect to the Earth equator. The mars parking orbit is 
assumed as 300 km circular and inclined at 75 degree with 
respect to the Mars equator. The high inclinations are cho- 
sen to satisfy the inclination constraint (Eq. (6)). The SOI 
duration of Earth and Mars are fixed as 3 and 2 days 
respectively after some trial runs. The trial runs indicate 
that lesser durations result in infeasible scenarios and lar- 
ger durations do not significantly improve the design in 
terms of achieved target parameters. Therefore, a suitable 
intermediate value is chosen to represent the SOI. The 
JPL ephemeris DE 405 is used to obtain the position and 
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velocity vectors of the planets. Vasile et al. (2005) and Olds 
and Kluever (2007) used evolutionary optimization tech- 
niques to generate preliminary trajectory design for 
Earth-Mars round trip and multiple gravity assist transfers. 
They used the concept of conventional patched conic tech- 
nique for the transfer trajectory design. In this paper, 
which deals with a direct transfer, a typical minimum 
energy opportunity is obtained using a simple grid search. 
A search over a range of departure epochs and flight dura- 
tions resulted in a minimum energy opportunity of 12 May 
2018 00:00:00 UTC with a flight duration 204 days and this 
opportunity is used for comparison and analysis purposes. 

The departure excess velocity vector with respect to 
earth equator and the arrival excess velocity vector with 
respect to mars equator are given in Table 2. These are 
obtained from the conventional and iterative patched conic 
methods. The excess velocity vectors for conventional and 
V-infinity tuned patched conic methods are same. While 
these methods generate only one excess velocity vector at 
the departure and arrival ends, the iterative patched conic 
method generates four distinct excess velocity vectors (cf. 
Table 2). 

The hyperbolic orbit characteristics obtained from the 
conventional and V-infinity tuned patched conic methods 


orbits in both the methods. Also, the design options 1-1 
and 2-1 have the same arrival hyperbolic orbits while the 
departure conditions are entirely different. Both these 
methods do not identify the distinct design options due 
to the absence of an iterative process that ensures continu- 
ity at the patch points. 

The hyperbolic orbit characteristics obtained from the 
proposed iterative method is given in Table 4. Note that 
the design options are distinctly different. The difference 
between the departure RAAN (Q,,,,) of design options 1- 
1 and 1-2 is 0.06 degree and the corresponding departure 
AoP (@.,) is 0.07 degree. These small differences in the 
departure angles, viz. RAAN and AoP, result in entirely 
different arrival geometries. This indicates the high sensitiv- 
ity of the arrival geometry to the departure angles. Similar 
trend is observed for the design options 2-1 and 2-2. Evi- 
dently, the proposed iterative patched conic method could 
identify the four distinct design options available for a 


Table 4 
Design options from the iterative patched conic method (i,.» = 75 deg, 
ix4 = 75 deg). 


é . . Parameters Option 1-1 Option 1-2 Option 2-1 Option 2-2 
are given in Table 3. The departure and arrival RAAN ; scons ; ood ; Tr ve ; 
of both the conventional and V-infinity tuned patched oe ae qitsa “ ‘ie 4 ian. 11 ot. 
conic designs are same. This is because, in V-infinity tuned = g, (deg) 333.3889 333.4465 129.9454 129.934] 
patched conic method, the analytical tuning that targets for — wp (deg) 167.3782 167.3057 64.4571 64.5547 
the V-infinity/excess velocity vector changes only the in- oo (km) —4980.0 —4981.3 —4975.7 —4977.1 
plane parameters of the hyperbolic orbit. Note that the  ¢~4 1.74240 1.74220 1.74304 1.74284 

: : : Qoo4 (deg) 68.0878 242.9041 68.0925 243.0652 

design options 1-1 and 1-2 have the same departure condi- 

: ; : : ; ‘ Mood (deg) 115.1783 314.9085 115.4581 314.5994 
tions but they result in entirely different arrival hyperbolic 
Table 2 
Departure and arrival excess velocity vectors for the minimum energy opportunity. 
Parameters Conventional/V-infinity Iterative patched conic method 

tuned patched conic method Option 1-1 Option 1-2 Option 2-1 Option 2-2 
Voop (km/s) 2.7891 2.7826 2.7839 2.7779 2.7791 
Xop (deg) 321.42 321.65 321.68 321.53 321.55 
Ovop (deg) —36.85 —37.21 —37.29 —36.84 —36.92 
VoocA (km/s) 2.9621 2.9602 2.9598 2.9614 2.9610 
Xo4 (deg) 245.66 245.51 245.49 245.60 245.56 
Ova (deg) 9.25 9.50 9.54 9.22 9.26 
Table 3 
Design options from patched conic methods (i, = 75 deg, ix4 = 75 deg). 
Parameters Conventional patched conic designs V-infinity tuned patched conic designs 
Option 1-1 Option 1-2 Option 2-1 Option 2-2 Option 1-1 Option 1-2 Option 2-1 Option 2-2 

Asxop (km) —51239.9 —51239.9 —51239.9 —51239.9 —58640.5 —58640.5 —58640.5 —58640.5 
CoD 1.13033 1.13033 1.13033 1.13033 1.11388 1.11388 1.11388 1.11388 
Qoop (deg) 333.0131 333.0131 129.8392 129.8392 333.0131 333.0131 129.8392 129.8392 
Moop (deg) 169.3999 169.3999 64.5845 64.5845 167.8129 167.8129 64.5845 64.5845 
Aso (km) —4881.1 —4881.1 —4881.1 —4881.1 —4973.4 —4973.4 —4973.4 —4973.4 
Coo 1.75744 1.75744 1.75744 1.75744 1.74339 1.74339 1.74339 1.74339 
Qooa (deg) 68.1673 243.1616 68.1673 243.1616 68.1673 243.1616 68.1673 243.1616 
Mooa (deg) 115.0951 314.2668 115.0951 314.2668 115.4118 314.5835 115.4118 314.5835 
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Table 5 


Comparison of analytical designs of option 1-1 with numerical design (i,.5 = 75 deg, i.4 = 75 deg). 


Parameters Conventional patched V-infinity tuned patched 
conic design conic design 

dsop (km) —51239.9 —58640.5 

Cx0D 1.13033 1.11388 

Qexop (deg) 333.0131 333.0131 

Mop (deg) 169.3999 167.8129 

dooa (km) —4881.1 —4973.4 

C04 1.75744 1.74339 

Qo (deg) 68.1673 68.1673 

Moa (deg) 115.0951 115.4118 


Iterative patched 


Iterative pseudostate Numerical design (realistic 


conic design design force model) 
—58965.7 —58625.7 —58104.1 
1.11325 1.11391 1.11493 
333.3889 333.6273 333.5550 
167.3782 167.5207 167.4887 
—4980.0 —4977.4 —4977.6 
1.74240 1.74270 1.74258 
68.0878 68.0867 68.0839 
115.1780 115.1685 115.1738 


direct interplanetary transfer by capturing the small differ- 
ences in the departure angles. 

Table 5 provides the different analytical designs includ- 
ing the iterative pseudostate design along with the numer- 
ical design for option 1-1. The numerical design is 
obtained by numerical search and numerical integration 
of equations of motion under a realistic force model. The 
realistic force model used in this paper includes: (1) Earth, 
Earth J5, Sun and Moon for 3 days during geocentric 
phase, (ii) Sun for 199 days during heliocentric phase and 
(iii) Mars and Sun for 2 days during the arrival planetocen- 
tric phase. Another force model, Patched Conic Force 
Model (PCFM) is used to evaluate the patched conic 
designs. The PCFM includes: (i) Earth for 3 days, (1i) 
Sun for 199 days and (iii) Mars for 2 days. 

The analytical patched conic designs (option 1-1) are 
numerically integrated under the patched conic force model 
and the achieved target parameters are given in Table 6. 
The desired closest approach altitude (CAA) at the arrival 
is 300km and the desired time of arrival periapsis is 2 
December 2018 00:00:00 UTC. The numerical integration 
of the conventional patched conic design resulted in a 
CAA of 2,736,169 km. The achieved time of periapsis is 
also largely deviated from the expected value. These large 
deviations occur because the SOI of the planets are not 
quantified and the hyperbolic orbit characteristics are not 
tuned to achieve the desired excess velocity vectors at the 
SOI. The achieved CAA from the V-infinity tuned patched 
conic design and iterative patched conic design are 
365,000 km and 205 km respectively. The deviations in 
the time of periapsis from these designs are about 3 months 
and 10 min respectively. There is a large improvement in 
the achieved CAA and the time of arrival periapsis with 
the iterative patched conic design. This is attributed to 
the fact that the iterative process connects the departure 


Table 6 
Evaluation of patched conic designs under force model PCFM. 


and arrival phases with the heliocentric phase and ensures 
continuity. Also, after each iteration on the patch point, 
the hyperbolic orbit characteristics are tuned to obtain 
the updated excess velocity vector at the SOI. 

The iterative patched conic design is also evaluated 
under the realistic force model and given in Table 7. The 
other analytical patched conic designs are not evaluated 
under the realistic force model because they do not gener- 
ate distinct designs. For comparison purposes, the evalua- 
tion of iterative pseudostate design under the realistic force 
model is given in Table 7. The proposed design achieves a 
CAA of about 464,872 km and the deviation in the time of 
periapsis is about 50h. The iterative pseudostate design 
achieves a CAA of about 265,194km and the deviation 
in the time of periapsis is about 44 h. The achievable accu- 
racies with the pseudostate design is better (but of the same 
order) because it includes the Sun’s effect within the planet 
SOI in the design process. Further, the evaluation of the 
performance of the iterative methods is carried out using 
two more parameters (cf. Table 7): (1) number of iterations 
for numerical refinement under the realistic force model 
when the analytical designs are used as the initial guess, 
and (ii) the trajectory correction maneuvers (TCM) 
required for the analytical designs to achieve the desired 
target conditions under the realistic force model. The ana- 
lytical design is numerically integrated and the target con- 
ditions are achieved by applying TCM at Earth’s SOI. The 
number of iterations and TCM requirement of the iterative 
pseudostate design is less compared to the iterative patched 
conic design. However, the iterative patched conic method 
derives its merit from its simple concepts and capability to 
identify the distinct design options. It is evident that the 
iterative patched conic design can also be a good initial 
guess for numerical refinement. The proposed method pro- 
vides additional information on the arrival angles such as 


Parameters 


Achieved CAA (km) 


Time of periapsis (UTC) 
DD MM YYYY HH:MM:SS 


Conventional patched conic design 
V-infinity tuned patched conic design 
Iterative patched conic design 205 


2,736,169 
365,000 


10 Mar 2019 00:00:00 
10 Mar 2019 00:00:00 
2 Dec 2018 00:09:15 
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Table 7 

Evaluation of iterative analytical designs under realistic force model, 

Parameters Achieved CAA (km) Time of periapsis (UTC) No. of iterations for numerical refinement TCM (m/s) 
DD MM YYYY HH:MM:SS 

Iterative patched conic design 464,872 4 Dec 2018 02:54:52 19 (40.7 s) 17 

Iterative pseudostate design 265,194 3 Dec 2018 20:01:45 14 (30s) 13 


RAAN which is essential to generate a particular arrival 
geometry in the numerical refinement process. In the 
absence of the additional information, the numerical refine- 
ment process may converge to any one of the possible arri- 
val geometries. 

To assess the computational efficiency of the proposed 
method, the computation time taken to generate designs 
for 2190 cases using different analytical methods are esti- 
mated in an Intel Core 15-3230 CPU 2.60 GHz processor. 
The designs are generated for a range of departure epochs. 
The departure epoch is varied over two years with a step 
size of 0.5 days, resulting in 2190 cases. The flight duration 
is fixed as 204 days. The time taken to generate these 
designs using the conventional patched conic and the V- 
infinity tuned patched conic methods are 204 and 220 ms 
respectively. The iterative patched conic method takes 
about 592 ms for a threshold value (difference between 
two successive patch points) of 1km. The proposed 
method generates improved designs at the small cost of 
computation time. Also, there is only a marginal increase 
in computation time for reduced threshold values. The per- 
formance of the proposed iterative method in terms of (1) 
number of iterations taken by the method to converge 
and (ii) the achievable accuracy under numerical integra- 
tion with the proposed design under PCFM force model, 
are shown in Table 8. A threshold value of 1 km for the 
patch points is used for design analysis purposes based 
on Table 8. 


5.2. Design analysis 


The FORTRAN code is used as a design analysis tool to 
analyse various mission scenarios. 


5.2.1. For different closest approach/APO periapsis altitudes 

The distinct transfer trajectory design options for differ- 
ent APO periapsis altitudes (200-1000 km) are analyzed 
using the iterative patched conic method. The APO apoap- 
sis altitude is fixed at 1000 km and the APO inclination is 
75 degree relative to Mars equator. The departure angles, 


Table 8 
Computation and precision accuracy of iterative patched conic method. 


Threshold value (km) Iterations Achieved CAA (km) 
1000 5 —114.565 

10 7 194.638 

1 10 205.263 

0.1 11 205.438 

0.01 12 205.461 


viz., RAAN and AoP of the design options are shown in 
Fig. 4. The variation in departure angles with APO periap- 
sis altitudes are marginal. The difference in departure 
RAAN between the options 1-1 and 1-2 is about 0.06 
degree for different APO altitudes. The difference in depar- 
ture RAAN between options 2-1 and 2-2 is even less, about 
0.012 degree. This shows that even a slight variation in the 
departure RAAN can result in entirely different arrival 
geometry. The difference in departure AoP between 
(options 1-1 and 1-2) and (options 2-1 and 2-2) are nearly 
same, about 0.07 degree. It is clear that the arrival condi- 
tions are extremely sensitive to the departure angles. 

The total velocity impulse (TPI+ POI) of the four 
design options for different APO periapsis altitudes are 
given in Fig. 5. There is no variation in the TPI impulse 
for the whole range of APO periapsis altitudes while the 
variation in POI impulse is about 100 m/s (cf. Fig. 5a). 
The total velocity impulses for the four design options 
for a particular APO altitude are almost same and so, the 
lines are not distinguishable for different APO altitudes 
as in Fig. 5a. Fig. 5b shows the variation for a shorter 
range of APO periapsis altitudes and it can be seen that 
the difference in total velocity impulse between different 
design options is only about 1 m/s. 


5.2.2. For different APO inclinations 

The distinct transfer trajectory design options are ana- 
lyzed for different APO inclinations. The departure date 
and flight duration are fixed as 12 May 2018 00:00:00 
UTC and 204 days. The declination of arrival V-infinity 
vector for this opportunity is 9.25 degree. Therefore, the 
APO inclination is varied from 10 degree to 170 degree. 
The variation of departure angles of the four distinct design 
options with different APO inclination are given in Fig. 6. 
The maximum difference in RAAN between options 1-1 
and 1-2 is about 0.07 degree and between options 2-1 and 
2-2 is about 0.03 degree. The maximum difference occurs 
near 90 degree inclination and the difference between the 
options is almost symmetric about this inclination. With 
a change of less than 0.05 degree in departure RAAN, 
the whole range of APO inclinations can be achieved. Sim- 
ilarly, the maximum difference in AoP between options 1-1 
and 1-2 is about 0.08 degree and between options 2-1 and 
2-2 is about 0.1 degree. Again, the maximums occur near 
90 degree inclination and the difference between options 
is almost symmetric about this inclination. The change in 
departure AoP required to achieve the whole range of 
APO inclinations is less than 0.05 degree. The sensitivity 
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Fig. 4. Departure angles of distinct designs for different APO periapsis altitudes (a) RAAN and (b) AoP. 
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Fig. 5. Total velocity impulse of four distinct design options for (a) a wide range of APO periapsis altitudes (b) for a short range of APO periapsis 
altitudes. 
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Fig. 6. Departure angles of distinct designs for different APO inclinations (a) RAAN and (b) AoP. 
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Fig. 7. Total velocity impulse of distinct designs for different APO 
inclinations. 


of the arrival geometry to the departure angles can be 
understood from this phenomena. 

Fig. 7 shows the variation of total velocity impulse of 
four distinct design options with different APO inclina- 
tions. The variation is symmetric about 90 degree inclina- 
tion. The maximum difference between options (1-1 and 
1-2) and options (2-1 and 2-2) are about 0.1 m/s and this 
occurs near 90 degree inclination. The variation in total 
velocity impulse for the entire range of APO inclinations, 
for all design options, is about 0.6-0.7 m/s only for a fixed 
opportunity. Note that the additional velocity impulse 
required to target any feasible APO inclination is negligibly 
small. 


6. Conclusions 


The proposed iterative patched conic method identifies 
the four distinct transfer trajectory design options available 
for a direct orbiter mission with a fixed departure epoch 
and flight duration. This method could capture the small 
differences in departure angles that result in entirely differ- 
ent arrival geometries. A design analysis tool is developed 
using the proposed method. The total velocity impulse 
required to achieve a particular arrival parking orbit 
(APO) periapsis altitude is almost same in all the four 
design options. The trans-planetary injection impulse 
required to achieve different APO periapsis altitudes is 
same. However, the additional total velocity impulse 
required to target any APO periapsis altitude varies to a 
maximum of about 100 m/s. For an opportunity, the addi- 
tional total velocity impulse required to target any feasible 
APO inclination is negligibly small. Also, with a change of 
less than 0.05 degree in departure Right Ascension of 
Ascending node (RAAN), the whole range of APO inclina- 
tions can be achieved. The change in departure Argument 
of Periapsis (AoP) required to achieve the whole range of 
APO inclinations is also less than 0.05 degree. The numer- 


ical propagation of the iterative patched conic design, 
under the patched conic force model, achieves a closest 
approach altitude of about 205 km, against the desired 
value of 300 km. Apart from providing a good initial guess, 
the proposed method provides additional information on 
arrival angles such as RAAN and AoP, which are required 
to generate distinct numerical designs. The proposed 
method retains the simplicity of the patched conic 
approach that involves only two-body models and identi- 
fies the four distinct design options available for a fixed 
opportunity. 
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